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1.  Introduction  and  Theory 


ALEGRA  is  a  finite-element  multiphysics  code  designed  for  modeling  shock 
hydrodynamics  and  coupled  electromagnetic  phenomena  including 
magnetohydrodynamics  (MHD).  This  multiphysics  capability  is  a  key  feature  of 
ALEGRA  and  the  result  of  many  years  of  multidisciplinary  effort.  Verification  and 
validation  (V&V)  of  ALEGRA  is  also  a  significant  undertaking.  Fortunately,  the 
proper  compartmentalization  is  ingrained  in  the  architecture  of  ALEGRA.  In  other 
words,  various  modules  of  ALEGRA  can  be  used  without  others  when  necessary. 
Therefore,  the  V&V  procedures  can  be  compartmentalized  as  well. 

We  have  2  goals  in  pursuing  this  project.  First,  we  explore  in  the  quasi-static 
approximation  the  evolution  of  the  magnetic  fields  inside  and  outside  an  inclusion 
and  the  parameters  for  which  the  quasi-static  approach  provides  for  self-consistent 
results.  Second,  we  explore  how  reliable  ALEGRA  is  in  its  static  limit,  specifically 
for  magnetic  diffusion.  By  the  static  limit  we  understand  the  stationary  states 
without  macroscopic  current.  We  choose  quite  a  general  class  of  3-D  solutions  for 
which  a  linear  isotropic  metallic  matrix  is  placed  inside  a  stationary  magnetic  field 
approaching  a  constant  value  at  infinity. 

Analysis  for  related  magnetic  diffusion  problems  appears  throughout  the  literature. 
Knoepfel1  considered  linear  and  nonlinear  magnetic  diffusion  for  simple 
geometries  and  nonpermeable  materials.  Rieben  and  White2  performed  simulations 
for  linear  transient  magnetic  diffusion  in  several  geometries,  including  annular  and 
spherical  shapes.  Woodson  and  Melcher3  analyzed  permeable  materials  for  a  slab 
geometry.  Brauer4  considered  slabs  and  cylinders  with  linear  and  nonlinear 
permeability  and  finite-element  modeling.  Here  we  consider  linear  permeability 
with  an  ellipsoidal  geometry  in  3-D.  This  expands  upon  the  previous  work 
described  in  Refs.  5  and  6,  which  examine  2-D  magnetic  diffusion  into  an  elliptic 
cylinder.  The  3-D  ellipsoidal  situation  is  considerably  more  challenging  because  of 
the  size  of  the  simulations  required  to  resolve  the  inclusion. 

1.1  MHD  Master  System 

The  analysis  of  quasi-statics  is  based  on  the  following  reduced  Maxwell  system: 

ldB1  ...  4? r  . 

=  ---jp  z‘l«VjHk  =  (1) 
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These  bulk  partial  differential  equations  (PDEs)  should  be  augmented  with  1)  the 
constitutive  equations  ]l  =  alJEj  (Ohm's  law),  2)  the  constitutive  equation  Bl  = 

Bl(Hk),  3)  boundary  conditions  =  0,  [/Z1]^  =  0,  4)  the  conditions  at 

infinity,  and  5)  appropriate  conditions  at  infinity  as  well  as  with  corresponding 
initial  conditions. 

Here,  zl>k  is  the  covariant  Levi-Civita  skew-symmetric  tensor;  £),  Ht,  and  Bt  are 
the  electric  and  magnetic  field  and  magnetic  induction,  respectively;  Jt  is  the 
electric  current  density  of  free  charges,  c  is  the  speed  of  light  in  vacuum,  and  is 

electrical  conductivity.  In  the  boundary  conditions,  n*  and  T*  are  the  normal  and 
tangent  vectors  to  the  discontinuity  boundaries.  The  ALEGRA  code  uses  the  vector 
potential  At.  The  vectors  At  and  Ht  are  interconnected  by  the  covariant  differential 
relation//'  =  zlikVjAk. 


1.2  Exact  Solution  for  an  Ellipsoidal  Inclusion 

There  are  few  exact  2-D  and  3-D  solutions  on  the  MHD  master  system.  For  the 
static  equilibrium  configuration,  a  closed-form  solution  can  be  obtained  for  an 
ellipsoidal  inclusion  in  an  infinite  isotropic  matrix,  in  particular,  in  vacuum,  which 
allows  generalization  for  anisotropic  and  nonlinear  media.7-9  Its  origin  traces  back 
to  the  pioneering  papers  by  Eshelby. 10,11  This  solution  is  described  in  the  following 
paragraphs  and  used  in  our  project  for  verification  purposes. 

Consider  an  ellipsoid  with  the  semi-axes  al5  a2,  and  a3  coinciding  with  the 
Cartesian  axes  z1,  z2,  and  z3.  We  assume  that  the  ellipsoidal  domain  is  filled  with 
a  linear  isotropic  substance  with  magnetic  permeability  /i.  We  then  assume  that  the 
ellipsoid  is  immersed  in  the  unbounded  space  in  which  there  is  a  uniform  magnetic 
field  Hio. 

If  there  is  an  ellipsoidal  inclusion,  the  otherwise  uniform  field  Hl  =  Hl°  will 
change.  The  changes  are  particularly  strong  inside  the  ellipse  and  in  its  vicinity.  At 
infinity,  the  newly  generated  field  Hl  approaches  its  original  value  Hl°. 

For  the  time-independent  fields  and  in  the  absence  of  macroscopic  currents  Jl,  the 
system  (Eq.  1)  splits  into  2:  one  for  the  electrostatic  field  and  the  other  for  the 
magnetic  field 


zijkVjHk  =  0  ,  (2) 

with  these  normal  and  tangential  boundary  conditions  at  the  surface  of  the 
inclusion: 


=  0, 
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lnl]+_n  =  o. 


(3) 


The  initial  conditions  become  unnecessary  in  this  case.  Let  Hl°  =  Bl°  be  the 
uniform  magnetic  field  in  vacuum  occupying  the  whole  space.  Consider  a  magneto¬ 
sensitive  body  a)  immersed  in  a  uniform  field.  When  the  body  is  introduced,  the 
magnetic  field  Hl°  and  induction  Bl°  change,  both  inside  and  outside  the  inclusion. 
If  the  inclusion  is  finite  in  size,  then  the  outside  field  Hl  =  Bl  in  the  vacuum  only 
asymptotically  approaches  the  uniform  fields  Hl°  =  Bl°  at  zl  oo  : 

W|*Hoo  =  fifzHoo  -  Hio  =  Bio,  (4) 

where  zl  are  the  spatial  coordinates.  The  disturbance  (Hl  —  Hl°)  will  be  generated 
because  of  the  dipole  magnetization  Ml  appearing  inside  the  inclusion,  but  whose 
influence  will  be  sensed  both  inside  and  outside  the  inclusion. 

This  problem  for  an  ellipsoidal  inclusion  in  unbounded  space  has  been  analyzed  by 
many  outstanding  mathematicians  and  physicists.  First,  the  problem  was 
considered  in  the  contexts  of  gravitation  and  cosmology.  Later  on,  the  problem  and 
its  solutions  found  multiple  applications  in  many  other  disciplines,  including 
electromagnetism,  and  it  is  in  this  context  that  we  use  it  here. 

We  seek  solutions  assuming  that  Hi  is  a  potential  field,  that  is. 

Hi  =  -Vtf  ,  (5) 

where  77  (z)  is  a  scalar  potential  field  that  should  not  be  confused  with  the  scalar 
potential  of  the  electrostatic  field,  or  with  the  vector  potential  A1  of  the  magnetic 
field.  When  dealing  with  an  ellipsoid  immersed  in  the  uniform  field  Ht  =  H°  at 
infinity,  we  will  be  looking  for  a  solution  77  _  (z)  inside  the  ellipsoid  in  the  following 
form: 


?7—  (z)  =  - Kiz 1 ,  (6) 

which  automatically  satisfies  the  PDE  in  Eq.  2  and  implies  that  Ht  =  Kt  inside  the 
ellipsoid. 

For  the  solution  outside  the  ellipsoidal  inclusion,  we  will  be  seeking  an  expression 
of  the  form 


77  +(z)=5iVi0+-Hi°zi, 


(7) 


where  0  is  the  Newtonian  potential  of  the  ellipsoid,  given  by  the  relationship 


dco* 

\z-z*\ 


(8) 
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and  it  satisfies  the  PDEs: 


VjVl0  =  —An  Inside  the  ellipsoid. 

VjV‘0  =  0  Outside  the  ellipsoid.  (9) 

Solutions  of  this  form  automatically  satisfy  the  PDE  in  Eq.  2  and  the  condition  at 
infinity 

//;(z)  -»  Hi°  at  |z|  ->  oo  ,  (10) 

and  it  is  implied  that  the  following  expression  holds  outside  the  ellipsoid: 

Ht  =  Hi°  -  G(VjQ  ,  (11) 

where  G?  is  a  uniform  tensor  to  be  determined. 

As  is  well-known,12  the  Newtonian  potential  within  the  ellipsoid  is  described  by 
the  quadratic  form 

0_(z)  =  C-^;zV,  (12) 

where  C  is  a  constant  and  YLj  is  a  symmetric  tensor.  Ytj  depends  only  on  the 
geometry  of  the  ellipsoid — we  call  it  the  “geometric  tensor”.  It  is  described  in 
Section  1.3.  Thus,  to  find  an  exact  solution  for  Ht,  we  need  to  find  the  constant  C 
and  the  2  unknown  vectors  Kt  (interior)  and  S;  (exterior). 

The  potential  0+(z)  outside  the  ellipsoid  is  much  more  complex.  Fortunately,  if  the 
potential  inside  the  inclusion  is  known  (i.e.,  the  constants  C  and  E^-  are  known), 
then  we  can  uniquely  recover  the  potential  outside  the  inclusion.  Using  the 
boundary  conditions  in  Eq.  3,  we  can  find  all  of  the  constants. 

Using  the  potential-based  definition  of  Ht  in  Eq.  5,  we  can  rewrite  the  surface 
tangential  and  normal  boundary  conditions  in  Eq.  3  as  follows: 

T]_=r]+  (tangential  component)  (13) 

and 

=  Vj?7+nl  (normal  component) .  (14) 

Then,  with  the  help  of  Eqs.  6  and  7,  we  get 

K.zi  =  — 5fcVfe0+  +  Hfz1 
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(15) 


and 


/r_A)rF  =  (— 5feVjVk0+  +  Hf0)  nl  .  (16) 

First  derivatives  of  the  potential  0  remain  continuous  across  the  ellipsoid  surface. 
Therefore,  we  can  rewrite  the  surface  tangential  boundary  condition  in  Eq.  15  as 

Ki  -  YkiSk  =  Hf  .  (17) 

Similarly,  we  can  rewrite  the  pointwise  surface  normal  boundary  condition  in  Eq. 
16  as  the  following  algebraic  condition: 

H-Ki+Sk(4n8ik-Yik)  =  Hi°,  (18) 

where  Sik  is  the  Kronecker  delta  tensor.  Using  Eq.17,  we  can  rewrite  Eq.  18  as 

(li.-l)Ki  +  4nSi  =  0.  (19) 

To  summarize,  Eqs.  17  and  19  comprise  a  system  of  2  linear  vector  equations  with 
2  unknown  vectors  and  Sj.  Eliminating  St  from  the  system  yields 

(ski  +  Yki^^)Kk  =  Hi°.  (20) 

After  solving  Eq.  20  with  respect  to  Kk,  we  can  find  St  from  the  equation 

/r_  —  1 

St=—K,  (21) 

Since  we  have  thus  found  the  vectors  Kt  and  5;,  the  determination  of  the  exact 
solution  is  complete,  and  all  that  remains  is  to  work  out  closed-form  solutions  for 
Hi,  which  can  be  used  in  verification  studies  for  ALEGRA. 


1.3  Geometric  Tensor 


To  build  closed-form  solutions,  first  the  geometric  tensor  Ytj  must  be  determined. 
The  main  components  of  Ytj  for  an  ellipsoid  with  semi-axes  equal  to  at  are  given 
by  the  formulae 


Ya  —  2na1a2a3  J* 
o 


dq 


(af  +  q)yj  (a?  +  q)(a|  +  q)(a|  +  q) 


(22) 


and  the  constant  C  appearing  in  Eq.  12  is  given  by 
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oo 


C  =  a1a2a3  I  ,  .  (23) 

l  V(ai  +  <?)(ai  +  <7)01  +  0) 

Here,  we  simplify  the  problem  by  limiting  our  interest  to  an  ellipsoidal  inclusion 
that  has  2  planes  of  symmetry — that  is,  it  has  2  axes  equal  in  length,  in  which  case 
it  is  a  spheroid.  In  the  case  of  a  spheroid  with  semi-axes  ax  =  a2  =  M  and  a3  =  TV, 
we  get 


Ylt  =  Y2  2  =  2nM2N  - .  and 

l  (M2  +  q)ViV2  +  q 


Y33  =  4n-  2Ytl 


(24) 


The  integral  in  Eq.  22  can  be  evaluated  using  elliptic  integrals.  For  a  spheroidal 
shape,  the  elliptic  integrals  reduce  to  the  more  elementary  trigonometric  functions. 
Defining  a  spheroid  aspect  ratio  e  =  N/M,  we  can  generate  expressions  for  oblate 
(e  <  1)  and  elongated  (e  >  1)  spheroids.  For  the  oblate  (saucer-like)  case  with  the 
axis  oriented  in  the  3-direction,  we  get 


Tn  =  Y22  =  2n  -  -  *2)3/2  (arccos  e  -  ejl^e2) 


and 


>33  =  4tt  -  2Ylt,  e<l.  (25) 

For  the  prolate  or  “elongated”  (cigar-like)  case  with  the  axis  oriented  in  the 
1 -direction,  we  get 


>22  =  >33  =  27 r  J  3/2  ( eje 2  -  1  -  arccos  e)  and 

Y1±  =  4? t-  2 Y22,  e  >  1  .  (26) 

The  geometric  tensor  can  be  written  more  generally  for  a  spheroid  as 

>fci  =  Yaxiallkh  +  (27T  —  “Taxiai)  (^fci  —  hh)  >  (2V) 


where  the  unit  vector  is  aligned  with  the  spheroid  axis,  and  Yaxial  is  given  by  T33 
in  Eq.  25  for  the  oblate  case,  and  by  Y1±  in  Eq.  26  for  the  prolate  case.  The  off- 
diagonal  components  of  the  geometric  tensor  are  zero  by  the  assumed  alignment  of 
the  spheroid  to  the  1,  2,  3  axes. 
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1.4  Closed-Form  Solution 


To  work  out  a  closed-form  solution  for  the  field  in  the  ellipsoid  interior,  we  first 
solve  Eq.  20  for  Kt.  Component-wise,  this  yields 

t  =  1 

k= 1 

1(1  i1  +  Yll^C~)  =  Hl° '  (28) 

i  =  2 

Yj{S«2Kt  +  Yj-^±K*)  =  H2° 

k= 1 

K2  (l  +  Y22  =  H2°  ■  (29) 

i  =  2 


fc  =  l 


K3  1  +  Y, 


33  ' 


47T 

~  1 
47T 


)=^3° 


)  =  #3°  ■ 


(30) 


Simplifying  these  expressions,  we  have 


Kl  =  Hio(l  + 


—  1 
47T 


(31) 


where  the  subscript  II  notation  indicates  diagonal  components  rather  than 
summation  over  repeated  indices. 

Having  found  Kt,  we  can  write  an  expression  for  the  field  in  the  ellipsoid  interior. 
To  obtain  the  exterior  solution,  we  must  also  substitute  Kt  into  Eq.  21,  as  follows: 


Si 


(32) 
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However,  the  form  of  the  potential  0+ outside  the  ellipsoid  is  not  analytically 
simple.  Thus,  even  having  found  St,  it  will  not  be  straightforward  to  comprehend 
the  behavior  of  Ht  outside  the  ellipsoid.  At  this  point,  we  proceed  by  limiting  the 
scope  of  our  study  to  the  ellipsoid  interior. 

The  ellipsoid  interior  has  a  relatively  simple  solution,  which  is  found  by  combining 
Eqs.  5  and  6  and  to  yield  Hi  =  A),  or,  from  Eq.  31, 


(33) 


Extending  this  to  the  magnetic  induction  and  retaining  the  assumed  orientation  of 
the  ellipsoid  axis  along  the  1 -direction,  we  have 


(34) 


Together  with  the  definitions  of  the  geometric  tensor  Yn  in  Eqs.  25  and  26  for  the 
oblate  and  prolate  spheroids,  we  have  a  simple  closed-form  solution  for  the 
magnetic  induction  interior  to  the  spheroid,  which  lends  itself  well  to  verification 
of  computational  electromagnetics  simulations.  This  represents  the  magnetic 
induction  remaining  in  the  material  after  magnetic  diffusion  has  saturated  and 
currents  have  decayed  to  zero,  leaving  the  steady-state  “equilibrium”  configuration 
of  the  field.  This  field  is  uniform  within  the  ellipsoid  and  depends  only  on  the 
geometry  of  the  ellipsoid,  its  relative  permeability,  and  the  imposed  field  at  infinity. 

2.  Numerical  Model  and  Simulations 


Having  now  arrived  at  a  transparent  closed-form  solution  for  the  magnetic  field  at 
equilibrium  inside  an  ellipsoidal  inclusion,  we  can  use  it  as  a  verification  tool  for 
numerical  modeling.  Transient  magnetic  diffusion  calculations  can  be  performed 
using  the  ALEGRA-MHD  code,13  and  the  accuracy  of  these  calculations  can  be 
probed. 

2.1  The  ALEGRA-MHD  Code 

The  “transient  magnetics”  module  of  the  ALEGRA-MHD  code  (henceforward 
ALEGRA)  computes  solutions  to  the  reduced  Maxwell  system  of  Eq.  1  in 
quasi-static  fashion.  It  is  assumed  the  medium  is  stationary,  with  an  electrical 
conductivity  a  and  a  magnetic  permeability  /i  that  are  fixed  for  each  material.  The 
system  is  recast  in  terms  of  the  vector  potential  At  and  transformed  to  SI  units,  and 

Approved  for  public  release;  distribution  is  unlimited. 


8 


appropriate  constitutive  relationships  are  incorporated.  These  include  Ohm’s  law 
Jl  =  aEl  and  a  simple  linear  relationship  between  the  magnetic  field  and  the 
magnetic  induction,  Bl  =  /J.H1.  An  implicit  linear  solver  is  used  with  a  finite- 
element  discretization  to  evolve  the  solution  forward  in  time. 

ALEGRA  is  equipped  to  handle  a  much  broader  class  of  problems,  including 
high- strain-rate  deformation,  mechanical  and  electromagnetic  forces,  and  Ohmic 
heating.  These  are  encompassed  within  ALEGRA’ s  broader  MHD  capability.  In 
these  systems,  body-fitted  meshes  and  purely  Lagrangian  approaches  are  usually 
impractical  because  of  large  strains  that  are  encountered.  Therefore,  ALEGRA 
includes  a  2-step  Lagrange-remap  formulation  that  permits  material  motion  across 
the  mesh,  and  the  presence  of  multiple  materials  in  a  single  element.  For  the  present 
work,  only  the  transient  magnetics  module  is  considered — that  is,  ALEGRA’ s 
capability  to  simulate  the  evolution  of  magnetic  fields  in  a  domain  containing 
stationary  materials  of  differing  conductivity  and  permeability.  But  connections  to 
the  broader  class  of  problems  are  retained  by  using  a  regular  Eulerian  mesh,  as 
would  be  used  in  those  situations,  with  material  interfaces  only  resolved  in 
volumetrically  mixed  multimaterial  elements.  This  facilitates  investigation  into 
whether  transient  magnetics  problems  can  be  accurately  represented  in  this  way. 
Considering  only  Eulerian  meshes,  however,  poses  a  significant  obstacle,  since 
material  interfaces  are  not  resolved  explicitly.  Here  we  examine  the  impact  of  the 
Eulerian  approximation. 

2.2  Verification  Strategy 

The  equilibrium  solution  for  the  ellipsoid  problem  described  previously  is  obtained 
in  ALEGRA  in  quasi-static  fashion  by  a  series  of  implicitly  integrated  timesteps 
capturing  the  time  evolution.  After  some  time,  the  transient  magnetic  diffusion 
process  saturates,  leading  to  a  distribution  of  the  magnetic  induction  that  is  no 
longer  varying  in  time.  At  this  time,  the  distribution  of  the  magnetic  induction 
exterior  to  the  inclusion  remains  spatially  nonuniform,  but  interior  to  the  ellipsoid, 
the  analysis  outlined  in  Section  1.4  shows  that  the  magnetic  induction  must  be 
spatially  uniform.  A  useful  verification  test  then  is  to  compare  the  configuration  of 
the  magnetic  induction  interior  to  the  ellipsoid  computed  by  ALEGRA  to  that 
predicted  by  the  analysis. 

To  carry  out  this  verification  test,  a  series  of  simulations  is  designed  for  ALEGRA 
in  3-D.  Dimensions  and  material  properties  including  conductivity  and 
permeability  are  chosen  for  the  inclusion  and  the  surrounding  region,  and  an 
imposed  magnetic  field  Hl°  is  defined.  Discretizations  are  then  chosen  to  allow  a 
convergence  analysis  by  spatial  refinement  of  the  mesh.  The  simulations  are  set  to 


Approved  for  public  release;  distribution  is  unlimited. 

9 


run  sufficiently  late  in  time  to  reach  the  equilibrium  condition,  and  comparison  to 
the  analytic  result  is  made  at  this  time. 

2.3  Sample  Simulation  Results 

Results  from  a  sample  3-D  ALEGRA  simulation  are  shown  in  Fig.  1.  Here,  an 
ellipsoid  with  a  major/minor  radius  of  1.8/0.56  cm  is  modeled  for  an  imposed 
magnetic  induction  of  1.0  T  oriented  in  the  y-direction  (vertical  in  these  images). 
The  ellipsoid  axis  lies  along  the  z-direction,  and  60  mesh  elements  span  the 
--dimension  of  the  ellipsoid.  It  has  a  conductivity  of  107  S/m  and  a  permeability 
exceeding  the  permeability  of  free  space  /i0  by  a  factor  of  3  (i.e.,  a  relative 
permeability  n_  =  3).  The  exterior  region  has  a  conductivity  of  10"6  S/m  and 
permeability  equal  to  /t0  (i.e.,  a  relative  permeability  of  1).  The  magnetic  induction 
Rj  evolves  over  time  until  a  steady  state  is  reached  and  negligible  current  persists. 
At  that  time,  the  magnetic  induction  inside  the  ellipsoid  is  uniform,  as  predicted  by 
the  preceding  analysis. 


F^s&j^OCOldr  d _ PsSvdOCOlOr  G - — - 

Van  BE_Y  t  -  £.0*04  hc  Var;  3E_¥  t  =  4.  Ce  >04  sec 


Fig.  1  Computed  solutions  for  the  vertical  magnetic  induction  ( By )  from  ALEGRA  during 

and  after  equilibration.  The  analytic  solution  predicts  1.5788  T  inside  the  ellipsoid. 

Evaluating  Eq.  34  for  this  geometry  and  imposed  magnetic  field  predicts  a  vertical 
magnetic  induction  By  inside  the  ellipsoid  of  1.5788  T.  In  the  simulation,  the  mean 
equilibrium  interior  field  is  approximately  By  =  1.596  T.  Comparison  with  the 
analytic  solution  reveals  that  the  approximation  used  to  insert  the  ellipsoidal  shape 
into  the  domain  must  be  improved  for  precise  comparison.  Increasing  the  number 
of  analytical  points  in  the  elliptical  cross  section  when  computing  initial 
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element-wise  volume  fractions  from  90  to  720  improves  the  computed  result.  The 
lateral  extent  of  the  simulation  mesh  also  affects  the  accuracy  of  the  result  and  must 
be  a  factor  of  3  to  4  larger  than  the  extent  of  the  ellipsoid  in  any  direction — only  a 
subregion  of  the  domain  is  shown  in  Fig.  1 .  Of  course,  the  mesh  resolution  affects 
the  accuracy  as  well,  and  this  is  explored  in  the  following  section. 

The  3-D  simulations  shown  in  Fig.  1  are  analogous  to  the  2-D  ALEGRA 
simulations  described  in  Refs.  5  and  6,  in  that  the  aspect  ratio  of  the  2-D  ellipse  has 
been  used  here  for  a  3-D  ellipsoid,  and  the  orientation  of  the  magnetic  field, 
dimensions,  and  material  properties  are  all  the  same.  However,  the  physical  time 
required  to  reach  a  steady-state  condition  is  much  smaller  in  this  3-D  ellipsoid  case 
than  in  the  2-D  ellipse  case.  For  the  ellipse  or  elliptical  cylinder  case,  roughly  2  ms 
is  required  for  the  rate  of  electromagnetic  energy  deposition  in  the  mesh  by  the 
imposed  magnetic  field  to  drop  by  3  orders  of  magnitude  from  its  initial  rate.  In  the 
3-D  ellipsoid  case,  only  about  0.5  ms  is  required  for  the  same  drop.  Magnetic 
diffusion  proceeds  to  a  steady  state  much  more  quickly  for  ellipsoid  than  for  the 
elliptical  cylinder  of  the  same  cross  section. 

2.4  Verification  Analysis 

To  provide  a  rigorous  verification  test  for  ALEGRA  in  3-D,  a  more  challenging 
geometry  than  that  shown  in  Fig.  1  is  used  here.  The  ellipsoid  is  elongated  to  an 
aspect  ratio  e  =  10,  and  an  oblique  imposed  magnetic  field  is  used,  given  by  Hl°  = 

(/r0V 3)  (x  +  y  +  z)  A/m.  This  geometry  is  shown  schematically  in  Fig.  2. 


The  ellipsoid  dimensions  (semi-axes)  for  the  verification  problem  are  a  =  3.16  cm 
in  the  x-direction  and  b  =  0.316  cm  in  the  y-  and  --directions.  Computational  mesh 
dimensions  are  w  =  24  cm  in  the  axial  (x)  direction  and  h  =  15  cm  in  the  transverse 
( y  and  z)  directions.  Based  on  the  outcome  of  the  verification  study  in  Refs.  5  and 
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6,  only  a  regular  rectangular  mesh  is  used  here,  with  no  body-fitting.  The  mesh  is 
composed  of  unit-aspect-ratio  (cubic)  hexahedral  elements  in  the  region  of  the 
ellipsoid,  with  gradually  larger  rectangular  elements  stretching  out  to  the 
boundaries. 

Convergence  of  the  steady-state  solution  at  t  =  1 .2  ms  under  spatial  refinement  is 
studied  by  simulating  this  geometry  in  ALEGRA  with  meshes  ranging  from  N=  20 
to  N  =  640  elements  spanning  the  axial  length  of  the  ellipsoid.  The  exact  solution 
predicted  by  Eq.  34  for  the  ellipsoid-interior  magnetic  induction  in  this  system  is 
Bx  =  1.6645184  T  and  By  =  Bz  =  0.8748995  T.  Error  estimates  for  these 
simulations  are  obtained  by  taking  an  L2  norm  over  all  values  of  Bx,  By,  and  Bz 
interior  to  the  ellipsoid  with  respect  to  this  uniform  exact  solution.  A  fractional  error 
norm  is  then  the  ratio  of  this  L2  error  norm  to  the  analytical  value  of  the  magnetic 
induction.  This  norm  can  be  measured  either  with  or  without  the  transitional  layer 
of  mixed-material  elements  that  span  the  void-material  interface  at  the  ellipsoid 
surface. 

Measurements  of  the  fractional  L2  error  in  the  solution  for  Bx  with  and  without  the 
mixed  elements  are  shown  in  Fig.  3.  We  see  that  the  error  is  dominated  by  the 
interfacial  region  where  mixed-material  elements  are  present.  When  the  interfacial 
region  is  included  (all  nonexterior  elements,  or  elements  where  ellipsoid  material 
is  present  in  any  amount),  the  error  magnitudes  are  larger  by  more  than  an  order  of 
magnitude,  and  the  convergence  rate  drops  from  roughly  1  to  roughly  0.5.  Data  for 
N  =  20  without  mixed-material  elements  are  not  included,  because  the  N  =  20  mesh 
is  so  coarse  that  all  elements  interior  to  the  ellipsoid  are  mixed  elements. 


Elements  spanning  spheroid  major  axis 

Fig.  3  Convergence  behavior  with  and  without  the  layer  of  mixed  elements  on  the  ellipsoid 
surface 
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Since  this  initial  measurement  demonstrates  that  the  mixed-material  elements  in  the 
simulation  carry  large  errors,  the  analysis  is  repeated  using  successively  smaller 
interior  subregions  of  the  ellipsoid,  with  lateral  and  axial  dimensions  75%  and  15% 
of  the  actual  inclusion.  The  convergence  trends  for  Bx  and  By  in  these  subregions 
are  shown  in  Fig.  4.  The  convergence  trend  for  Bz  is  nearly  identical  to  that  for  By, 
so  it  is  not  shown  here. 


Elements  spanning  spheroid  major  axis 


Elements  spanning  spheroid  major  axis 


Fig.  4  Convergence  trends  for  (a)  the  x-component  and  ( b )  the  j-component  of  magnetic 
induction,  for  various  subregions  of  the  spheroid  inclusion 

The  solution  computed  by  ALEGRA  converges  monotonically  at  approximately  first 
order  to  the  predicted  analytic  solution.  The  convergence  rate  and  error  magnitudes 
improve  as  the  ellipsoid  sampling  subregion  shrinks,  and  fewer  elements  near  the 
interface  are  included  in  the  error  sum.  The  error  magnitudes  are  significantly  smaller 
for  Bx  than  for  By  and  Bz,  because  of  the  major/minor  axis  asymmetry  that  has  been 
built  into  this  test.  The  convergence  rates  are  shown  in  Table  1. 

Table  1  Convergence  rates  for  spheroid  verification  study  for  various  subregions. 


Case 

Rate  for  Bx 

Rate  for  By 

Rate  for  Bz 

All  nonexterior  elements 

0.525 

0.365 

0.365 

No  mixed  elements 

0.755 

0.675 

0.675 

Central  75% 

0.864 

0.759 

0.759 

Central  15% 

1.013 

0.961 

0.961 

2.5  Origin  of  Diminished  Convergence  Rates 

Since  verification  is  done  here  using  magnetic  induction  B  (of  most  interest  to 
ALEGRA  users),  rather  than  the  native  vector  potential  A,  we  expect  only 
first-order  convergence  of  the  solutions.  This  is  because  the  magnetic  induction  is 
not  the  native  finite-element  solution  variable,  and  ALEGRA  must  use  certain 
approximations  to  obtain  B.  Thus,  first-order  convergence  with  respect  to  spatial 
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mesh  refinement  is  what  we  expect,  and  that  is  what  we  observe  in  the  tests,  so  long 
as  the  region  near  the  material  interface  is  excluded. 

Further  investigation  into  the  origins  of  the  larger  error  magnitudes  and  diminished 
convergence  rates  near  the  interface  shows  that  anomalous  spatial  oscillations 
persist  in  the  solution  even  at  steady  state.  These  are  shown  in  Fig.  5  for  the  3  finest 
mesh  resolutions,  where  spatial  distributions  of  the  fractional  difference  between 
the  computed  and  analytic  values  of  Bx  are  plotted  on  a  slice  through  the  ellipsoid 
centroid  normal  to  the  z-direction.  These  oscillations  also  appear  in  the  2-D 
solutions  discussed  in  Ref.  5  for  the  regular-mesh  cases. 

The  oscillations  are  unphysical  and  arise  because  in  the  mixed-material  elements, 
the  natural  interface  conditions  on  B  and  H  (see  Eq.  3)  are  not  enforced  explicitly, 
since  the  interface  is  actually  treated  as  a  volumetrically  mixed  zone  at  the  level  of 
the  Eulerian  spatial  discretization.  Further,  when  an  ALEGRA  simulation  includes 
materials  with  distinct  magnetic  permeabilities,  the  element-level  homogenization 
scheme  actually  averages  the  reluctivity,  not  the  permeability.  Such  prominent 
oscillations  have  not  been  observed  in  magnetic  diffusion  simulations  where  the 
magnetic  permeability  is  /t0  everywhere.  Therefore,  these  results  suggest  that  the 
option  to  form  an  average  element  permeability  might  be  useful. 


N=  160 


N  =  320 


N  =  640 


Fractional  error  in  Bx 

-1  .Oe-02  -2.5e-03  5.0e-03  1 ,2e-02  2.0e-02 


Fig.  5  Spatial  distribution  of  error  in  ellipsoid-interior  magnetic  induction  at  steady  state, 
arranged  by  increasing  number  of  N  elements  spanning  the  axial  length  of  the  ellipsoid 
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2.6  Computational  Cost 


The  ALEGRA  simulations  conducted  here  vary  enormously  in  computational  cost. 
At  N  =  40,  the  mesh  has  16,800  elements,  and  the  simulation  runs  to  completion  in 
a  few  minutes  on  8  processors.  At  N=  640,  the  mesh  has  56  million  elements  (1.7 
billion  edges)  and  runs  to  completion  in  about  3  days  on  4,096  processors.  This  is 
very  near  the  upper  limit  of  the  computational  domain  size  that  can  be  used  in 
ALEGRA. 

3.  Conclusions 


This  work  shows  how  the  equilibrium/steady-state  condition  of  the  magnetic  field 
in  an  ellipsoidal  inclusion  can  be  obtained  analytically.  The  analysis  results  in  a 
simple,  elegant  closed-form  expression  that  describes  a  uniform  magnetic  field  in 
the  ellipsoid  interior  exactly.  The  expression  is  easily  evaluated  within  verification 
testing  techniques,  making  it  well  suited  for  use  in  verification  of  electromagnetics 
modeling  tools  like  ALEGRA. 

Convergence  testing  using  this  analytic  solution  shows  that  the  equilibrium 
magnetized  state  can  be  reached  by  transient  means  via  computation  with 
ALEGRA.  The  solution  in  the  interior  core  of  the  ellipsoid  converges  to  the  exact 
solution  at  first  order,  as  expected,  for  a  very  large  range  of  spatial  mesh  sizes 
spanning  the  very  coarsest  to  the  very  finest  meshes  that  are  possible  using 
ALEGRA. 

The  verification  work  reveals  2  significant  issues  in  the  computation  of  magnetic 
diffusion  for  magnetically  permeable  materials  on  Eulerian  (nonbody-fitted) 
meshes.  Those  issues  are  1)  that  the  natural  boundary  conditions  on  B  and  H  are 
not  enforced  explicitly  and  2)  that  element-level  homogenization  of  the  elements 
on  the  ellipsoid  interface  using  reluctivity  can  be  inaccurate.  These  issues  result  in 
noticeable  local  perturbations  from  the  exact  uniform  solution,  which  corrupt  the 
ellipsoid  interior  field  slightly  and  cause  deterioration  of  the  convergence  rate  when 
these  elements  are  included  in  the  error  metric. 

Verification  testing  is  applied  here  by  developing  an  analytic  solution  and  using 
convergence  analysis.  The  test  results  can  be  used  to  advance  the  algorithms  in 
ALEGRA  and  codes  like  it.  It  is  anticipated  that  future  work  extending  those 
algorithms  can  make  use  of  the  analytic  solution  and  results. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


2-D 

2-dimensional 

3-D 

3-dimensional 

MHD 

magnetohydrodynamics 

PDE 

partial  differential  equation 

V&V 

verification  and  validation 

Approved  for  public  release;  distribution  is  unlimited. 

17 


1  DEFENSE  TECHNICAL 
(PDF)  INFORMATION  CTR 

DTIC  OCA 

2  DIRECTOR 
(PDF)  USARL 

RDRL  CIO  L 

IMAL  HRA  MAIL  &  RECORDS 
MGMT 

1  GOVT  PRINT  G  OFC 
(PDF)  AMALHOTRA 

3  USARL 
(PDF)  RDRL  WMP  A 

J  CAZAMIAS 
RDRL  WMP  C 
M  GRINFELD 
RDRL  WMP  D 
R  DONEY 

2  SANDIA  NTL  LAB 
(PDF)  J  NIEDERHAUS 

C  SIEFERT 


Approved  for  public  release;  distribution  is  unlimited 


